Statistical inversion, also called Bayesian statistics, is a method for reconstructing an image from measurement data. Typically the image is a representation of the values of a parameter, like attenuation coefficient for a particular kind of radiation, in a large number of points of a three-dimensional target region, and the measurement data comes from a two-dimensional imaging detector. As an example, we may consider the tomographic X-ray imaging of a part of a human or animal body, although it must be remembered that 3D imaging and statistical inversion can be equally applied also to a large number of other purposes.
FIG. 1 illustrates a basic principle of computerised tomography with sparse angle data in two dimensions. An object 102 is irradiated from a limited number of directions with a radiation source, which in FIG. 1 is shown in two exemplary positions 106 and 107. A spatially sensitive detector, shown correspondingly in two positions 104 and 105, gives spatially bound attenuation information for each exposure. In FIG. 1 the circular regions 103 inside the object 102 are assumed to attenuate radiation more than the bulk of the object. Also, the region 101 outside the object is also assumed to have a different attenuation than the bulk of the object. The readings from the detector and the corresponding irradiation angles are taken into a computer 108, which uses a mathematical reconstruction method to calculate what sort of configuration of attenuation regions is located inside the object 101 in order to give just these particular attenuation profiles. The result of the calculation, which essentially represents an attenuation coefficient map of a two-dimensional slice of the object 102, is shown on a display 109. In order to gather additional information for the measurement, external measurement means like a stereoscopic camera 110, a mechanical measurement arm 111 or a laser scanner 112 can be used to scan the location and form of the outer boundary of the object 102.
Bayesian statistics is known to produce very good images for limited angle tomography reconstructions. Advances in mathematical methods, explained for example in the Finnish patent number FI 116324 and the Finnish patent application number FI 20125327 (which at the time of writing this description is not yet available to the public) have advocated Bayesian statistics as a practical method for 3D image reconstruction.
When formulating a forward model for tomography, one needs to parameterize the medium in some way. FIG. 2 shows an exemplary way of doing this, by dividing the modelled target region into equally spaced pixels (in the two-dimensional case) orvoxels (in the three-dimensional case, illustrated in FIG. 2). The forward theory of the propagation of the ray along the path of each ray passing from a radiation source to the detector can then be modelled using one of many possible interpolation methods. An example of a ray is shown as 201 in FIG. 2, and the detector is illustrated schematically as 202 in FIG. 2. The grid has an exaggeratedly coarse resolution in FIG. 2 for demonstration purposes only. In practice the grid may contain thousands or even millions of voxels: at the time of writing this description a three-dimensional grid in a medical application could have e.g. 512×512×512 or even 2000×2000×2000 voxels.
The core of Bayesian statistics is the so-called likelihood function p(m|x). Ignoring possible calibration parameters, the likelihood function p(m|x) of a tomographic model is completely defined by the set of model parameters xi, which define the two- or three-dimensional map of attenuation coefficients (or other characteristic descriptors) of the medium. The subscript i refers here in general to the way in which the model parameters are identified. For example in FIG. 2 a three-dimensional indexing scheme is utilised, with the parameters x111 to x144 of the 16 front voxels explicitly illustrated. Pixels of the detector are similarly indexed, so that in FIG. 2 we have a two-dimensional detector of M×N pixels and an exemplary two-dimensional indexing scheme mmn of detector pixels, where the indices mε[1,M] and nε[1,N]. At the time of writing this description a large 2-dimensional X-ray detector may have 1000×1000 pixels.
In addition to Bayesian inversion, other approaches are known to the reconstruction of 3D images from 2D measurements. One such approach is known as the simultaneous algebraic reconstruction technique, or SART. A detailed description of a SART-based approach can be found for example in Hengyong Yu and Ge Wang: “SART-Type Image Reconstruction from a Limited Number of Projections with the Sparsity Constraint”, International Journal of Biomedical Imaging, Volume 2010 (2010), Article ID 934847, 9 pages, doi:10.1155/2010/934847.
Independent of whether the mathematical approach is based on Bayesian inversion or SART, a set of detection results can be expressed in vector form asm=Ax+ε  (1)where m is a vector of detection results, x is a vector of unknown variables (i.e. the parameter values of interest at all voxels of the grid), and ε is an error term. The so-called theory matrix A contains elements that indicate, what is the contribution of each voxel in the grid to a detection result obtained at a particular pixel of the detector. Since there may be e.g. 10002 pixels in the detector (and consequently 10002 components in the vector m), and e.g. 20003 voxels in the grid, the theory matrix A becomes very large: in this exemplary case it has 1·106 rows and 8·109 columns, which makes 8·1016 elements in total for one projection only.
Allocating two bytes for storing the value of each element in such a theory matrix would require 16 petabytes (16·1015 B) of memory, which at the time of writing this description is beyond the capability of most computer systems or would at least require an exceptionally large and expensive memory. The theory matrix is relatively sparse, and only its nonzero elements need to be stored, but the required memory space is still easily at least hundreds of gigabytes. The size of detectors and the resolution required of tomographic models may be expected to increase in the future, and many imaging cases require a large number of projections, which means a corresponding increase in memory requirements. It seems that the required amount of available memory as well as the speedy access to stored data may become major bottlenecks in applying Bayesian statistics to real-life imaging tasks.